#ifndef LMMBERMUDANSWAPTIONPRICER_HPP
#define LMMBERMUDANSWAPTIONPRICER_HPP
#include "LMMPricer.hpp"
#include "LMMPathPool.hpp"
#include "leastSquareSolver.hpp"

namespace TermStructure
{
	template<typename TVarianceStructure,template<typename, typename,typename> class TMultiDimensionGaussian, typename T = double,typename TSize = long>
	class LMMBermudanSwaptionPricer : public LMMPricer <TVarianceStructure, TMultiDimensionGaussian, T, TSize>
	{
	protected:
		LMMPathPool<TVarianceStructure, TMultiDimensionGaussian, T, TSize> _pathpool;
	public:
		T PriceBermudanSwaption(Swap<T,TSize> &contract, TSize NumberOfPath = 10000)
		{
			T result = 0;
			this->_pathpool.SetLMMAdvancer(this->UnderlyingAdvancer);
			this->UnderlyingAdvancer.AnchorIndex = contract.SwapEndIndex;
			this->_pathpool.NumberOfPath = NumberOfPath;
			this->_pathpool.CurrentTimeIndex = 0;
			this->_pathpool.EndIndex = contract.SwapEndIndex - 1;
			this->_pathpool.Generate();
			SimpleMatrix<T,TSize> cvx(NumberOfPath, 6);
			SimpleMatrix<T,TSize> cvy(NumberOfPath,2);
			SimpleMatrix<T,TSize> coef(6,1);
			Swap<T,TSize> anotherContract(this->UnderlyingAdvancer.GetTimeSequence(), contract.SwapEndIndex, 0);
			anotherContract.Strike = contract.Strike;
			anotherContract.SwapEnterIndex = contract.SwapEndIndex  - 1;
			for(TSize i = 0; i < NumberOfPath; i++)
			{
				cvy(i,0) = anotherContract.Price(this->_pathpool(i, contract.SwapEndIndex - 1), contract.SwapEndIndex - 1);
				if(cvy(i,0)<0) cvy(i,0) = 0;
				else cvy(i,0) /= this->_pathpool(i,contract.SwapEndIndex - 1).DiscountFactorAt(contract.SwapEndIndex, contract.SwapEndIndex - 1);

			}
			
			for(TSize i = contract.SwapEndIndex - 2; i > contract.SwapEnterIndex; i--)
			{
				//std::ofstream out("tester.csv");
				anotherContract.SwapEnterIndex = i;
				for(TSize pi = 0; pi < NumberOfPath; pi++)
				{
					cvy(pi,1) = anotherContract.Price(this->_pathpool(pi,i),i);
					cvy(pi,0) *= this->_pathpool(pi,i).DiscountFactorAt(contract.SwapEndIndex, i);
					T x = _pathpool(pi,i).SwapRate(i+1, contract.SwapEndIndex);
					T y = _pathpool(pi,i)[i+1];
					cvx(pi,0) = 1;
					cvx(pi,1) = x;
					cvx(pi,2) = y;
					cvx(pi,3) = x*y;
					cvx(pi,4) = x*x;
					cvx(pi,5) = y*y;
					//out<<x<<","<<y<<","<<cvy(pi,0)<<std::endl;
				}
				LeastSquareSolver<T,TSize>::Fit(cvx,cvy,coef);
				//outputTo(cout,coef);
				for(TSize pi = 0; pi < NumberOfPath; pi++)
				{
					T temp = 0;
					for(TSize l = 0; l < 6; l++)
					{
						temp += cvx(pi,l)*coef(l,0);
					}
					if(temp < 0) temp = 0;
					cvy(pi,0) = temp > cvy(pi,1) ? temp: cvy(pi,1);
					cvy(pi,0)/= this->_pathpool(pi,i).DiscountFactorAt(contract.SwapEndIndex, i);
				}
				//out.close();
			}
			for(TSize pi = 0; pi < NumberOfPath; pi++)
			{
				result += cvy(pi,0)/NumberOfPath;
			}
			_pathpool(0,0).EstablishDiscountFactors();
			return result*_pathpool(0,0).DiscountFactorAt(contract.SwapEndIndex, 0);
		}
	};
}

#endif